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ABSTRACT 

Numerical calculations of the linear Rossby wave instability (RWI) in global three-dimensional 
(3D) disks are presented. The linearized fluid equations are solved for vertically stratified, radially 
structured disks with either a locally isothermal or polytropic equation of state, by decomposing 
the vertical dependence of the perturbed hydrodynamic quantities into Hermite and Gegenbauer 
polynomials, respectively. It is confirmed that the RWI operates in 3D. For perturbations with 
vertical dependence assumed above, there is little difference in growth rates between 3D and two- 
dimensional (2D) calculations. Comparison between 2D and 3D solutions of this type suggest the 
RWI is predominantly a 2D instability and that three-dimensional effects, such as vertical motion, 
to be interpreted as a perturbative consequence of the dominant 2D flow. The vertical flow around 
co-rotation, where vortex-formation is expected, is examined. In locally isothermal disks the expected 
vortex center remains in approximate vertical hydrostatic equilibrium. For polytropic disks the vortex 
center has positive vertical velocity, whose magnitude increases with decreasing polytropic index n. 



1. INTRODUCTION 

Theoretical modeling of protoplanetary disks lead to 
complex structures that are unlikely to be described 
by s mooth radial proflles (jTerquemI l2008t lArmitagg 
1201 Ih . However, radially structure d disks may develop 
the Rossby wa ve instability (RWI, iLovelace et al.|[T999l : 
iLi et al.l 120000 , w hich leads to v ortex-formation in the 
nonlinear regime ()Li et al.ll2001|) . Thus, the RWI may 
play a role in the evolution of protoplanetary disks. 

The disk RWI is a dynamical instability associated 
with the presence of extrema in the ratio of vorticity 
to surface density, or vortensitjQ. The instability results 
from wave coupling across such an extremum. Its physics 
is similar to the Papaloizou-Pri ngle i nstability (PPI, 
Paoaloizou & Prinalc 1984, 198^, [l987l : iGoldreich et"aLl 
1986t INaravan et al.l 119871 ) which operate in pressure- 
supported, thick tori. The RWI operates in thin, 
centrifugally-supported disks with non-power law rota- 
tion proflles, and is insensitive to radial boundary condi- 
tions. 

The relevance of the RWI in protoplanetary disks has 
been demonstrated in two situations. Variability in 
the efficiency of turbulent angular mo mentum transport 
by th e magneto-rotational instability (jBalbus fc Hawlev 
19911) can result in the existence of 'dead zones' (jGammie 
1996f). in which the turbulent viscosity is small. The ra- 
dial boundary between a dead zone and the actively ac- 
cretin g region is prone to the RWI (Varnicr e fc Taggen 
[200l ILvra eTall VEM . [20091: ICrespe ct al. 201l|), with 
observable consequences (jRegalv et al. 2012). In ad- 
dition to hydrodynamic angular momentum transport, 
the RWI may also assist planet formation forma- 
tion by co ncentrating solids into anti-cyclonic vortices 
(jBarge fc S ommeria 1995). 

Another origin of the RWI in protoplanetary disks. 
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^ This quantity is modified by a factor involving the disk entropy, 
if the latter is not constant. 



which motivated this study, is disk-planet interaction 
(jGoldreich fc Tremainel[l979l. [l980il . A sufficiently mas- 
sive planet leads to gap opening (jLin fc Papaloizoul 
119861 ). while low mass protoplanets may op en gaps pro- 
vided the disk viscosity is sufficiently small (jMuto et al.l 
120101 : iDong et 311120111 ) . Vortensity jumps across planet- 
induced shocks lead to the necessary disk profile for the 
RWI to develo p around gap edges (iKoUer et al.l 120031 : 
ILi et all l2005t Ide Val-Borro et al.l 120071 ). Subsequent 
vortex-formatio n significantly affects disk-planet torques 
and migration COu et al."2007t ILi et all 120091 : lYu et al.l 
2010; Lin & Papaloizou 201(|T^ 

The above studies of the RWI have all employed the 
razor-thin disk approximation (but not e that the PPI 
was originally analyzed in 3D). lYu fc Lil ([2009) have ex- 
amined the RWI with a toroidal magnetic field in a 3D 
but unstratified disk. iMeheut et al.l ()2010[ ) first demon- 
strated the RWI in nonlinear hydrodynamic simulations 
of 3D stratified disks ( later with improved resolution 
in Mc heut et al]l2012a[ ). while lUmur had (|2010D analyzed 
the RWI in approximate 3D disk models. 

Recently, IMeheut et"all ()2012b[ ) calculated linear RWI 
modes in a three-dimensional, globally isothermal disk, 
which displayed vertical motion. In this paper, we com- 
pute linear RWI modes in three-dimensional disks across 
a range of parameter values, including different equations 
of state. Our focus here is on how these affect the vertical 
flow in the co-ro tation region, w here vortex-formation is 
known to occur (|Li et al.ll200lD . 

This paper is organized as follows. In f}2] we list the 
governing equations and describe our disk models. We 
derive the linearized fluid equations in fj3] and describe 
our numerical methods in [g] Results are presented in 
fj5] for locally isothermal disks and in fj6] for polytropic 
disks. In fjT] we briefly examine a linear mode qualita- 
tively different to those above, found in a disk model in- 
volving K^ < (taken from lMeheut et al.l l2010). where k 
is the epicycle frequency. We summarize and discuss our 
results in ^Hl including limitations of our calculations. 
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2. GOVERNING EQUATIONS, DISK MODELS AND 
ASSUMPTIONS 

We consider a three-dimensional, inviscid, non-self- 
gravitating disk orbiting a star of mass M* and adopt 
(r, (j), z) cylindrical polar co-ordinates centered on the 
star. The frame is non-rotating. The governing equa- 
tions are the 3D Euler equations: 



|+V.(p.) = 0, 

dv 1 

-— +v -^v = --VP - V<E>, 

at p 

P = P{r,p), 



(1) 

(2) 
(3) 



where p is the density, P is the pressure, v is the velocity 
field and $* is the gravitational potential due to the cen- 
tral star. Eq. [3] is an equation of state (EOS), specified 
later, such that the pressure may be calculated without 
an energy equation. 

We assume the disk is geometrically thin so that $* 
may be approximated as 



$*(r,z) = - 



GAf, 



GAU 



^j,2 _|_ ^2 



1 



2r2 / 



(4) 



This approximation is adopted so that the resulting equi- 
librium density field has a convenient functional form 
suitable for the application of orthogonal polynomials 
(see 35). This greatly simplifies the numerical prob- 
lem. Henceforth we use the approximate potential for 
self-consistency. 

The unperturbed disk is steady, axisymmetric with no 
meridional velocity {dt = d^j, = Vr = Vz =0). The disk 
is stratified with p = p{r, z) set by vertical hydrostatic 
balance. The azimuthal velocity is Vff, = rJl, where J7 is 
the angular speed. Vff, is set by radial balance between 
pressure, stellar gravity and centrifugal forces. Because 
the disk is thin, the angular velocity is close to Keplerian, 
i.e. 17 ~ f]fe = {GM^/r^f/'^. 

To introduce radial structure, we choose the unper- 
turbed surface density profile to be 



E(r) 



pdz 



= Sn 



To 



l + {A-l) exp 



{r - r^f 
2Ar2 



(5) 



(|Li et al.ll2000t ). Eq. [S] corresponds to a Gaussian surface 
density bump centered at r = ro, width Ar and ampli- 
tude A, on top of a background power-law profile with 
index —a. Since disk self- gravity is ignored, the surface 
density scale Sq is arbitrary. 

To specify the three-dimensional disk structure, we 
choose the EOS to be either locally isothermal or poly- 
tropic. These are described below. 

2.1. Locally isothermal disks 
For locally isothermal disks the pressure is calculated 



as 



P = cl{r)p, 



(6) 



where Cs{r) is the sound-speed given by c^ ~ Hflk and 
H{r) is the disk scale-height. The unperturbed density 
is 



Pir,z) = 



E(r) 



27riJ(r) 



cxp 



2H^{r) 



(7) 



In the numerical calculations we will choose H{r) — hr 
with h being a constant aspect-ratio, since this is a typ- 
ical model for protoplanetary diskto. The exponential 
decay means the gas density becomes negligible after a 
few scale heights. Thus the vertical domain can be taken 
to be z g [— oo, oo], even though we have made the thin- 
disk approximation. 

2.1.1. Approximate equilibrium 
For simplicity, we set the azimuthal velocity to 



rdP 
p dr 



a$. 



z=0 



dr 



(8) 



z=0 



Away from the midplane the deviation from exact ra- 
dial balance is proportional to h^ <^ 1 for a thin disk 
(jTanakaet al.|[2062[ ). We adopt Eq. [8] to allow us to 
apply standard solution methods. 

To test whether or not the precise form of Vt affect 
our results, we also considered setting p — > S in Eq. [8l 
which gives the velocity profile W0.2D for a razor-thin disk. 
For our fiducial calculation, growth rates differ by ~ 1% 
between adopting Eq. [8] or W0,2Dj and we observe the 
same flow structure. 

In fact, locally isothermal disks generally have differ- 
ential rotation in z, i.e. i7 = ^{r,z), unless the disk 
is also globally isothermal. It is therefore important to 
note that in assuming Eq. [S] , we have artificially sup- 
pressed baroclinic effects. We discuss some justification 
for this in i)8.4l and Appendix [X] Although the chosen 
basic state is not in exact equilibrium, setting Vt = Vl(r) 
greatly simplifies the linear equations as the only vertical 
dependence of the basic state is through the exponential 
factor in p. It allows us to address the specific question 
of whether or not vertical density stratification has any 
effect on the RWI, without the complication of baroclin ic 
instabilities (|Knobloch fc Spruitl[T986l : lUmurhanll2012[ ). 

2.2. Poly tropic disks 

In order to set up a more self-consistent basic state, 
that is, Q — ^{r) and a finite vertical domain, we also 
consider polytropic disks, for which 



P = Kp'+i, 



(9) 



where K is a constant and n is the polytropic index. 
Vertical hydrostatic equilibrium imply 



p(r, z) 



GM^H^jr) 
2K{l + n)r^ 



m{r) 



Po(r) 



m{r) 



(10) 



^ This choice also enables us to compare the locally isothermal 
disk with a polytropic disk with constant aspect-ratio. 
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Here, z = H is the disk surface where p{r, H) — 0. Thus, 
when discussing polytropic disks H is referred to as the 
disk thickness. 

The function H{r) and mid-plane density po(^) Skie cal- 
culated through 



E(r)=po(r)iJ(r)/„, 



(11) 



where /„ = /_j(l — x^)"dx, with pQ{r) related to H{r) 
by Eq. [10] and I](r) given by Eq. [S] We can therefore 
write 



H{r) = Ho 



S(r) 



(12) 



where Hq = H{ro) is the disk thickness at the bump ra- 
dius. We parametrize it by writing Hq = hro so that h 
is the aspect-ratio at tq. Note that a surface density en- 
hancement by a factor A corresponds to an enhancement 
of the disk thickness by a factor ^^/(^"+^). 

For a polytropic dis k the azimuthal velocity is s trictly 
independent of z (e.g. IPapaloizou fc Pringlelll984D . It is 
given by 






1- 



3IP 

'2^ 



H dH 

r dr 



(13) 



where the second equality follows from the approxima- 
tion for the stellar potential in a thin disk (Eq. U). 

Of course, given H{r) one can obtain the azimuthal 
velocity v^^e corresponding to the exact gravitational po- 
tential of a point mass. For our fiducial setup, the differ- 
ence in growth rate is < 4% between using w^^e and using 
V0 above, and we observe no difference in flow structure. 
However, we will use Vij, so that the equilibrium density 
and velocity fields are self-consistent and in exact balance 
with the same potential. 

3. LINEARIZED EQUATIONS 

In this section we derive the governing equation for 
small disturbances in the disk. As described above, the 
basic state is p = p{r, z) and v = (0, ri7, 0), with il — 
fl(r). The perturbed state is assumed to have the form 

/9 ^ /9 + Re[(5p(r, z) expi((Ti + mi/))], (14) 

P^ P + Re[SP{r, z) expi{<Tt + m(j))], (15) 

V ^ V + Kc[Sv{r,z)expi{(Tt + m(p)], (16) 

where a ^ a^ + ij is a complex frequency {an, 7 being 
real) and m is the azimuthal wave-number taken to be a 
positive integer. We will omit writing 'Re' below, with 
the understanding that physical solutions correspond to 
real parts of the complex perturbations. 

For the locally isothermal equation of state, the lin- 
earized momentum equations give 



D 



SVr = ;^ 



5v, = 



dr 



2mmV\ 



K^ dW amW 



D \2n dr 



dz ' 



(17) 
(18) 
(19) 



where W = Sp/p is the relative density perturbation, 
a = (7 + mn{r) is the shifted frequency, D = k^ — a^, 
and 



2 " / 4o2\ 



(20) 



is the square of the epicycle frequency. Corresponding 
equations for the polytropic disk are very similar, and 
are readily obtained by setting Cs to unity and replacing 
W ^ S = SP/p where S is the enthalpy perturbation. 

Inserting the perturbed velocity field into the linearized 
continuity equation 

1 (J I'm u 

iaSp + -^ {rpSvr) H pSv^ + —- {pSv^) = 0, (21) 

r or r az 

yields, for locally isothermal disks: 



r6p 



d f rpc^ dW\ 2mW d f c^p^ 



dr \ D 



dr J 



9 9 
m CgP 

rD 



W 



a dr 
rci d 



D 



dW\ 



(22) 



and for polytropic disks: 

d frpdS 
dr \ D dr 



rSp 



2mS d /pf7\ 
a d^ V"d J 



m^p 



^_^d_ f dS_ 
CT^ dz V dz 



(23) 



We remark that Eq. [22] is in fact valid for locally 
isothermal disks with any fixed sound-speed profile Cs{r), 
assuming the equilibrium azimuthal velocity is indepen- 
dent of 2; ( Appendix l^l. Also note that Eq. [23] is actually 
valid for any barotropic EOS , i.e. whenever P — P{p). 
The 3D problem is to solve Eq. [22]— [23l which will gen- 
erally describe disturbances depending on (r, 0, z) and 
motion in all three directions. 

3.1. Relation to the two-dimensional problem 

We define the 2D problem as solving Eq. [22]— [23] sub- 
ject to dz — 0. Denoting the corresponding solutions 
as W2D {r) , 5*20 (r) and inserting them into the governing 
equations yields, after vertical integration, 

r3^ ^A. ( r^cldW2D \ ^ 2mW2D d |^c|Sf7^| 



dr \ D 



dr 



D 



- 



rD 



SM^2D 



(24) 



for locally isothermal disks and 

_ d /rEdS'2D\ 2mS2D d /T.fl\ 
~d^ V"d dr ) d- d^ \D' ) 



rD 



S: 



2D 



(25) 



for polytropic disks, where ST, — J_ 6pdz is the sur- 
face density perturbation. Note that VF2D — ST/T, is the 
relative surface density perturbation, and 6*20 = ^H/E 
where Sli is the perturbation to the vertically integrated 
pressure (H = J^^Pdz). Solutions to Eq. [24]-[25] de- 
scribe disturbances which only depend on (r, 0) and there 
is no vertical motion. 
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As defined here, the 2D problem and 3D problem 
involves the same background disk, which is three- 
dimensional. However, the governing equation for lin- 
ear disturbances in razor-thin disks have the same form 
as Eq. B^l— B5] when the razor-thin disk has a locally 
isothermal or barotropic EOS in the form 11 = c'^{r)Y, or 
n — n(i;), respectively. 

3.2. Co-rotation singularity and the RWI 

Inspection of the 2D equations, Eq. [51}— HU reveal a 
potential singularity when <7{rc) — 0, where re is the 
co-rotation radius defined by 



aii + mfl{rc) — 0. 



(26) 



This co-rotation singularity can be rendered ineffective 
if rr also satisfies 



= for locally isothermal disks, (27) 

= for polytropic disks, (28) 

VIZ r^ 

where 



d 


\V J 


d 


fl\ 


dr 


\VJ r 



v 



2fiE 



(29) 



is the vortensity. T he qu antity ??/c^ can be seen as a gen- 
eralized vortensity ()Li e t al. 2000), but for convenience 
we will simply use 'vortensity' in the discussion below. 
Thus there can exist 2D neutral disturbances with co- 
rotation at a vortensity extremum, for which the 2D lin- 
ear operator is real and regular everywhere. 

Strictly speaking, co-rotation singularities only con- 
cern neutral disturbances (7 = 0). In practice we are 
interested in growing solutions (7 < 0) so such singular- 
ities do not arise in the numerical computation. Nev- 
ertheless, the discussion above is important because the 
growth rates we find are typically I7I ^ ^(rci). Further- 
more, association of r^ with a vortensity extremum forms 
the basis of the RWI. 

In studies employing razor-thin disks, the RWI has 
largest disturbance amplitude in the co-rotation region 
where |a^| <C k^. It can be shown that such modes 
can only be unst able if the re exists vortensity extrema 
in the disk (e.g. iLin fc Pap aloizou 2010). Indeed, the 
RWI is found to have wi th co-rotation radius r^ close 
to a vortensity minimum (iLovelace et al.lll999l : iLi et al.l 
l2000HLin fc Papaloizoull2011aD . 

It is precisely linear modes with the above properties 
which we wish to explore in 3D. However, we do not 
expect such modes to have significant 2-dependence in 
their relative density or enthalpy perturbation around co- 
rotation. From the linearized vertical equation of motion 
we see that 

IdX 

(7 OZ 

where X is W ot S depending on the EOS. Near co- 
rotation where \a\ is small, \dzX\ should be almost negli- 
gible. Otherwise, even small vertical gradients in density 
or enthalpy perturbation will cause significant vertical 
motion, and linearization becomes invalid. 



4. NUMERICAL PROCEDURE 

In principle one could attempt a numerical solution to 
the partial differential equations (PDE) above, for ex- 
ample by finite-differencing in the (r, z) plane. However, 
since one of our goals is to assess three-dimensional ef- 
fects, it is more useful to have a numerical scheme that 
automatically separates out the 2D problem from the full 
3D problem. 

We begin by making the co-ordinate transformation 



{f,z)^ 


{r,z/H), 


(30) 


d d\ 

dr' dz ~ 


[df "" H dz' Hdzj' 


(31) 



where ' denotes differentiation with respect to the argu- 
ment. In this co-ordinate system the background den- 
sity is separable, i.e. pir,z) = g{f)f{z), where / = 
exp (— z^/2) for locally isothermal disks and / = (1 — z^)" 
for polytropic disks. This motivates us to seek solutions 
of the form 



W^ = ^W^/(r)H/(z), 

1=0 

00 

S = Y,Si{f)C^{z), 



(32) 
(33) 



1=0 



where Hi is a Hermite polynomial of order / and C^" is 
a Gegenbauer polynomial of index A and order I. Note 
that radial and vertical variations are coupled because 
z = z{r) through H{r). 
These polynomials satisfy the orthogonality relations 



(34) 



nkiz)niiz)cxpi~z'/2)dz = V2T:ndki, 



(35) 

where Ski here is the Kronecker delta aii d F is the Gamma 
function (| Abramowitz fc StegunI [19651) . For polytropic 
disks, we choose the parameter A to be 



X = n . 

2 



(36) 



Consequently, for a polytropic index n = 1.5, C^ are 
the Chebyshev polynomials of the second kind, and for 

1/2 
n — 1, Cj are the Legendre polynomials. Eigenfunction 

expansions in z is a standard method to account for ver- 
tical dependence in disk pro blems (e.g. ^O kazaki fc Kate 
1985; Paoaloizou & Pringldl 119851 iTakeuchi fc Mivama 
1998; Tanaka et al...20oC ^ 

It is important to keep in mind that by assuming 
the above decompositions (Eq. [551— [55]) we restrict the 
type of perturbations to those satisfying certain physi- 
cal conditions implied by the orthogonal polynomials at 
the upper disk boundary. In the locally isothermal disk 
we require th e kinetic energy density to be bounded at 
large heights ([Takeuchi fc Mivamalll998l) . and for poly- 
tropic disks a regularity co ndition applies at z = ±1 
(iPapaloizo u fc Pringlelll985[ ). Such perturbations can be 
decomposed as above because the polynomials form a 
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complete set (jZhang fc Lail 120061 ) . On the other hand, 
the above specific decomposition cannot be applied if 
one considers other vertical boundary conditions (e.g., 
conditions imposed at other heights). 

After transforming the governing equations into (f, z) 
co-ordinates, we insert the anstaz Eq. l32]— i33l into Eq. 
\T2\ — Eq. [221 multiply by Hk, and C^ respectively, then 
integrate vertically. This procedure yields an equation of 
the form 



AiXi + BiXi^2 + CiXi+2 = 0, 



(37) 



where Xi is Wi or Si , and Ai, Bi, Ci are linear operators 
which only depend on r and a, but are different for the 
two EOS (see Appendixl^. For each I there is a separate 
equation with the operators Bi , Ci representing coupling 
with the / ± 2 modes. Note that Bi is set to zero when 
/ = 0, 1. 

We have now transformed the governing partial dif- 
ferential equation into an infinite set of coupled ordinary 
differential equations (ODE). In practice we truncate the 
solution at /max, i-e. Xi = for I > ?max- The decom- 
position has the advantage that for modes nearly inde- 
pendent of z, Imax can be small. In the simplest case of 
setting /max ~ 0, we only solve 

AoXo = 0, (38) 

which is the 2D problem. That is, if /max — then 
Wo = W21) and Sq = 5*20- 

4.1. Matrix methods 

We now proceed to a numerical solution to the linear 
problem. We discretize the linear operators and solutions 
on a grid which divides the radial range r G [?'i,ro] into 
Nr uniformly spaced points. The coupled set of ODEs 
then become a single matrix equation. This is denoted 
generically as 



Mx = 0, 



(39) 



where the square matrix M represents the discretized 
linear operator and the vector x is the discretized solu- 
tion. The size of the matrix and vector depends on /max- 
For example setting /max = 4, Eq. [SHlthen represents the 
discretized version of 

AoXo + C0X2 = 0, 

B2XQ + A2X2 + C2X4 = 0, 

B4X2 + AiXi = 0, 

for which M is a 3Nr x 3Nr matrix and a; is a vector of 
length 3Nr. 

The matrix problem, Eq. 1391 is a set of homogeneous 
linear equations. Non-trivial solutions exist if 



detM = 0. 



(40) 



The complex frequency a is required to be such that the 
matrix M{a) is singular. We have used two approaches 
to achieve this. The first is to consider the usual eigen- 
value problem: 



Starting with a trial a, standard matrix softwarqj may 
be used to find the eigenvalues p and associated eigen- 
vectors X. We then apply Newton- Raphson iteration to 
solve fmin/|t'max| = by Varying a, where fmin,max cor- 
responds to eigenvalues of smallest and largest absolute 
value found from Eq. 1411 

Another approach is to perform a singular value de- 
compositioi:Q (SVD) of M, so that 



M^Udmg{si,S2,...)V\ 



(42) 



M{a)x = vx. 



(41) 



where [/, V are unitary matrices (^ denotes Hermitian 
conjugate) and the real numbers Si > are the singular 
values of M. The columns of V are the right singular 
vectors of M . If min(si) = then Mxq = 0, where 
Xq is the right singular vector associated with min(si). 
We therefore use Newton-Raphson iteration to zero the 
quantity F = xIMxo/x'qXq by varying a. 

These methods give the same result. We always per- 
form the SVD for the final matrix M{a) in order to eval- 
uate R~^, where R = max(si)/min(si) is the condition 
number of M. Since R = oo for a singular matrix, we 
only accept solutions for which R^^ is zero at machine 
precision (typically R~^ < 10~^^). The matrix meth- 
ods ou tlined above was also used in iLin &: Papaloizoul 

(UnniS). 

4.2. Radial boundary conditions 

For simplicity we impose dXi/dr = at r = 
ri, Tq. The RWI is associated with internal structure 
away from boundaries. Consequently, it is insensi- 
tive to radial boundary conditions in razor-thin disks 
(|de Val-Borro et al.|[20Q7l: ILin fc Papaloizoull2011aD . We 
assume this still holds in 3D. For exa mple, approximate 
3D disk models developed bv lUmurha n (2008, 201_3), in 
which the inner/outer disk boundaries play no role, also 
support the RWI. 

As a check, additional calculations were performed 
with: drX = applied at boundaries (which introduces 
mode coupling), different r^, Tq and a numerical condi- 
tion where boundary derivatives are approximated by in- 
terior points. The last case is strictly a numerical proce- 
dure to generate a closed set of equations to solve. For 
the solutions of interest, these experiments gave results 
with no appreciable difference. 

4.3. Fiducial setup 

We work in units such that G — M* — 1. Our standard 
disk spans r G [?'z,?'o] = [0.4, 1.6] and has a surface den- 
sity profile with a = 0.5. The bump is located at ro = 1 
with width parameter /3 = 0.05. We use Nr = 512 grid 
points and first solve the 2D problem (/max = 0), then 
use the obtained eigenvalue to start the iteration for the 
3D problem, for which /max = 6. We only consider even 
/. 

In pl we will use the setup employed bv lMeheut et al.l 
()2010f rto examine a 3D RWI mode with k^ < Q aX the 
bump radius. This mode appears quite different to our 
standard setup with k^ > everywhere. 

3 We used LAPACK. 

■* We used LAPACK for a direct decomposition. We 

also performed the SVD with PRDPACK (available at 
http://soi. Stanford. edu/$\sim$rniunk/PROPACK/), which is 
an iterative method. I'hese gave the same results. 
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Fig. 1. — Background profile of tfic fiducial locally isothermal 
disk with A. = 1.25 and h = 0.07, in terms of the generalized 
vortensity 'rj/c^, scaled by its value at the bump. Unstable modes 
are associated with the minimum at rg. 

4.4. Results analysis 

The solution to the Unear problem gives the complex 
radial functions Xi , which can be used to reconstruct the 
complex amplitudes, e.g. Svz{r,z) by using Eq. 15^1— [551 
and Eq. I19[ but we are interested in physical (real) solu- 
tions. We will often visualize the solution for a specific 
m with two-dimensional plots. We explain below how 
these are obtained. 

The real perturbation is, e.g. Ke[Svz exp i(at + m(f>) , so 
the spatial dependence of a physical perturbation is 

Svz — i> Re[5wz] cos (to0) — Im[(5w2] sin {m<p), (43) 

and similarly for other variables. We focus on the so- 
lution near the vortex core, defined to be at (r, (p) = 
(ro,0o), where 



cos(m</)o)= Re[X{ro,Q)]/\X{ro,0)\, 
sin(m0o) - -Im[X(ro,0)]/|X(ro,0)|. 



(44) 



The magnitude of the (real) perturbation is arbitrary 
but its sign is fixed, e.g. X, now representing the real 
density or enthalpy perturbation, is positive at (r, </>, z) = 
(ro, 00 1 0). In practice the vortex core is near a maximum 
midplane over-density. 

We visualize results in the (r, z) plane by setting </> = 
00 in Eq- SSI Similarly, perturbations are visualized in 
the (r, (/)) plane at a chosen z, and in the (0, z) plane 
at r = ro with the azimuthal range set to € [0o ~ 
7r/2?7i, 00 + 7r/277i]. For convenience we also define VIq = 
ri(ro) and kq = ^(ro). 

5. RESULTS: LOCALLY ISOTHERMAL DISKS 

For locally isothermal disks we choose h — 0.07 and 
A = 1.25 as a fiducial case. Recall c^ ex 1/r, so that far 
from rg the generalized vortensity ij/c^ is fiat, and is a 
minimum at tq. The background disk is shown in Fig. [TJ 
Note that k^ > everywhere, and min(K^/r2|) ~ 0.59. 

Recall that for locally isothermal disks we assumed 
an approximate basic state ( §2.1.ip . The extent of in- 
exact radial balanc e in the background depends on h 
(jTanaka et al .1 [20021 ) . In a nonlinear simulation this may 
lead to radial motion. To keep this effect fixed in com- 
paring different linear calculations below, in this section 
we fix h. 

5.1. Solution example 



TABLI 

ElGENFREQUENCIES 
ISOTHERMAL DLSK 


] 1 

IN THE LOCALLY 
WITH h = 0.07 


"1 -(Tfl/(mQo) 


-lOS/f^o 


1 0.9960 (0.9960) 

2 0.9960 (0.9960) 

3 0.9961 (0.9960) 

4 0.9964 (0.9964) 

5 0.9972 (0.9971) 

6 0.9980 (0.9978) 


2.8038 (2.8044) 
4.8931 (4.8985) 
5.7205 (5.7365) 
5.1245 (5.1843) 
3.4557 (3.5720) 
1.8317 (1.9615) 



Note. — Values in brackets were ob- 
tained from the 2D problem. 

We solved the fiducial case for m S [1,6]. Table [T] com- 
pares the eigenfrequencies obtained from the 2D and 3D 
problems. Growth rates in 2D and 3D are very similar, 
so the instability is largely associated with Wq . We thus 
expect the RWI to grow in 3D disk s on the same time- 
scales as in the razor-thin disk^ fe.g. lLi et al.ll2000l ). The 
growth rate for the most unstable mode (m = 3) is only 
~ 0.06J7o but this corresponds to ~ 3 orbits at rg, so the 
instability operates on dynamical time-scales. 

Fig. [5] compares the radial functions Wi for the m — 3 
and m — 5 modes. In both cases. Wo dominates over 
Wi^o, implying that the relative density perturbation is 
nearly z-independent. For ?n = 3, Wq itself is dominated 
by the co- rotation region r ^ rg, but for m — 5 the 
amplitude in the oscillatory region is larger than that 
around ro. The Wi>2 modes are negligible, so three- 
dimensional effects are due to W2. Unlike Wq, in both 
cases IW2I has largest amplitudes in the wave-like regions 
towards the boundaries, and is smallest near rg. This is 
consistent with the absorption of waves with I > at 
co-rotation discussed in Liet al- (2003). 

It is well known that in the razor-thin disk, as m is 
increased the RWI becomes more w ave-like (as see n here 
for Wo) and is eventually quenched (jLi et al.ll2000[ ). This 
might contribute to the slightly smaller growth rates ob- 
tained in the 3D problem than in the 2D problem (Table 
[T|) since Wi^q are wave-like (in addition to wave absorp- 
tion at CO- rotation). However, this effect is unimportant 
because their amplitudes are much smaller than Wq. 

Although the observed stabilization effect increases 
with m, Wq loses its RWI character at high m. Thus, 
it can be said that the RWI, considered as a low m, 
radially confined non-axisymmetric disturbance, has a 
growth rate determined by the 2D problem. 

5.2. Flow in the (r, 0) plane 

Fig. [3| shows the perturbed velocity field in the (r, 0) 
plane for the m — S mode above, and for a case with 
A ^ 1.6 (growth rate ~ 0.15rio)- The figure shows that 
anti-cyclonic motion at an over-density is a robust fea- 
ture. This confirms that the unstable modes found here 
are indeed the analog of the RWI in razor-thin disks. 
The perturbed horizontal velocity {Svr,5vij,) has negligi- 
ble variation with respect to z. 

^ This statement assumes the 2D problem give similar growth 
rates to the equivalent razor-thin disk setup, which we have checked 
to be the case. 
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Fig. 2. — Radial eigenfunctions Wi for locally isothermal disk 
with h = 0.07. These are normalized by |Wo(r'o)|. The I > modes 
have also been magnified in order to compare its radial structure 
with Wo- 




Fig. 3. — Perturbed horizontal velocity field in the (?■,(/>) plane, 
for the locally isothermal disk with h = 0.07 and A = 1.25 (left) 
and with ji = 1.6 (right). The contours indicate relative density 
perturbation. The case with yl = 1.6 dis plays a double- peak in 
density perturbation, which is explained in lLi et al.l I I200CI ). 



5.3. Vertical flow 

We now examine vertical flow associated with the RWI. 
We focus on the co-rotation region since this is where 
relative density perturbations are largest and vortex- 
formation is expected. 

Fig. |4] shows the perturbed vertical velocity field in 




Fig. 4. — Vertical velocity field (contours) for the m = 3 mode 
in the locally isothermal disk with h = 0.07, in the (r, z) plane 
at azimuths 6 = m{ij> — ij>o) = 0.2tt (top), (middle) and — 0.27r 
(bottom). Arrows indicate the perturbed velocity field projected 
onto this plane. 

the (r, z) plane, at several azimuths. Since the largest 
contribution to Svz comes from the I = 2 term in the 
expansion for W (Eq. l32t , the magnitude of Svz increases 
linearly with z. 

Ahead and behind the vortex core, the flow just fol- 
lows the anti-cyclonic motion, with radial variations in 
Svz being negligible. At (p = (j>o there is also very little 
vertical motion for z < Q.bH, but there is upwards mo- 
tion at r = 0.9, 1.1, i.e. the edge of the vortex (see Fig. 
[3|). This can affect how dust particles are collected by 
RWI vortices. 

For comparison. Fig. [S] shows the vertical flow for the 
m = 2 mode. This flow is more two-dimensional than 
the fiducial case above. This is expected for decreasing 
m (se e, e.g. iPapaloizou fc Pringlelll985l : iGoldreich et al.l 
|1986[ ). It also appears qualitatively different (e.g. down- 
wards flow at r = 1.1 instead of upwards as see for 
TO = 3). We typically flnd locally isothermal disks to 
display a wider range of flow patterns around co-rotation 
than polytropic disks presented later, which show generic 
patterns. 

Finally, Fig. [B] shows the perturbed vertical velocity in 
the {4>tZ) plane at r = rg. Vertical motion is upwards 
ahead of an anti-cyclonic vortex and downwards behind 
it. The vertical velocity can be comparable to the per- 
turbed azimuthal velocity, so the perturbed flow is fully 
three-dimensional in this plane. However, the vortex cen- 
ter (ro, c^o) remains in vertical hydrostatic balance. This 
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Fig. 5. — Vertical velocity field (contours) for the m = 2 mode 
in the locally isothermal disk with h = 0.07. The slice is taken 
at the azimuth 4) = <j>o. This figure is to be compared with the 
middle plot in Fig. |4] Arrows indicate the perturbed velocity field 
projected onto this plane. 
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Fig. 6. — Vertical velocity field (contours) for the m = 3 mode 
in the locally isothermal disk with h = 0.07, in the (</», z) plane 
at the radius r = tq. Arrows indicate the perturbed velocity field 
projected onto this plane. 

is not the case for polytropic disks. 

5.4. Dependence of vertical flow on instability strength 

We now assess how the three-dimensionahty of the flow 
in the co-rotation region varies with instabihty strength. 
We examine ratio (|(5w2|)/(|(5tir|) , where (•) denotes aver- 
aging over r S [0.9, 1.1] and z £ [0, 2H] at fixed azimuth 
(f) = (f>o- In calculating this ratio, we ignore W/>2 be- 
cause the dominant contribution to Svr and Svz comes 
from Wo and W2 respectively. This ratio is large if there 
is significant vertical motion. 

Results are shown in Fig. [71 where the bump ampli- 
tude A is increased at fixed h = 0.07. Growth rates 
increase with A , which is expected from previous works 
(|Li et al.ll20Q0f ). but the flow actually becomes less three- 
dimensional with increasing instability strength. 

In the co-rotation region where a ^ ij, we expect from 
the linearized equation of motion that 



|(5dz 



H 



W2 

7 



n'. 



(45) 



\Svz\ scales with I/I7I, so that increasing growth rates 
contributes to decreasing \Svz\- Thus, the flow in the co- 
rotation region does not necessarily become more three- 
dimensional with increasing A. 

It is clear from Fig. [7] that three-dimensionality de- 
creases because of increasing I7I since (|W2|)/(|Wo|) 
varies weakly. We demonstrate this in Fig. |51 which 
shows that in the disk with A = 1.6 the flow is mainly 
horizontal. As in the fiducial case with A = 1.25, there 




Fig. 7. — Average magnitude of vertical velocity (solid), in the 
co-rotation region of the RWI in the locally isothermal disk, as a 
function of bump amplitude yl at fixed aspect-ratio h. Also shown 
are the normalized amplitude of W2 in this region (dashed) and 
the growth rates (dotted). 




Fig. 8. — Same as Fig. |4]but for a disk with A = 1.6. The slice 
is taken at (p = ipo. 

is little motion at r ^ ro. 

Fig. [71— [8] shows that in the locally isothermal disk, 
more unstable modes are also more two-dimensional (in 
the CO- rotation region). \W2\ remains a small fraction of 
I Wo I and \Svz\ is largely affected by |7|. 

However, I7I can be obtained by just solving the 2D 
problem. Thus, we could have anticipated the trend of 
\6vz\ in Fig. [71 based on only 2D calculations, with the 
assumption that changes in \W2\ are less significant than 
the increase in I7I . The above explicit calculation confirm 
this, suggesting we interpret the RWI as predominantly 
a 2D instability and that three-dimensional effects on the 
RWI are small (for low m). We further illustrate these 
points with polytropic disks below. 

6. RESULTS: POLYTROPIC DISKS 

Our fiducial polytropic disk has polytropic index n = 
1.5. In the absence of a bump, a surface density profile 
oc r^^/^ gives a constant aspect-ratio {H c>c r). The 
bump parameters are set to ^ = 1.4 and h = 0.14. Recall 
that for polytropic disks, H is the disk thickness and h 
is the aspect-ratio at tq. 

Although the surface density enhancement is relatively 
large, it corresponds to only ~ 9% enhancement of the 
disk thickness at tq. The background disk is shown in 
Fig. ini in terms of the vortensity profile. The fiducial 
disk has a global vortensity gradient (77 oc r~^ away from 
ro), but it is the local minimum that drives instability. 
The epicycle frequency is such that min(K^/rj^) = 0.47. 

6.1. Solution examples 




Fig. 9. — Background profile of tfie fiducial polytropic disk with 
n = 1.5, A = 1.4 and h = 0.14, in terms of the vortonsity. 
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Fig. 10. — Radial functions Si for the n = 1.5 fiducial polytropic 
disk. These are normalized by | So (ro ) | . The I > modes have also 
been magnified in order to compare its radial structure with So- 



TABLE 2 

ElGENFREQUENCIES IN THE n = 1.5 

POLYTROPIC DISK 



m 


-cru/irnQo) 


-lo27/f7o 


1 
2 
3 
4 
5 
6 


0.9930 (0.9930) 
0.9934 (0.9934) 
0.9941 (0.9941) 
0.9947 (0.9946) 
0.9952 (0.9945) 
0.9954 (0.9950) 


4.4900 (4.4907) 
8.2793 (8.2867) 
10.769 (10.793) 
11.594 (11.591) 
10.646 (10.861) 
8.0092 (8.5802) 



Note. — Values in brackets were ob- 
tained from the 2D problem. 

Eigenfrequencies for the fiducial case are shown in Ta- 
ble [2j The modes of interest are those with disturbance 
amplitudes largest near rg, which were found to corre- 
spond to m < A. These modes have effectively the same 
growth rates in 2D and 3D. This gives confidence that 
the RWI is an 2D instability. We will consider m = 3 
below in order to compare with locally isothermal disks. 
The m = 3 growth rate is only slightly smaller than 
the most unstable m = 4 mode. In co-rotation region, 
low m m odes are also insensitive to radial boundary 
conditions dLin fc Papaloizoull20lTah . 

Fig. [To] shows the Si functions for the m — 3 case. 
These are similar to the locally isothermal disk (Fig. [2]). 
We typically find the / > radial functions to have larger 
amplitudes (compared to ^ = 0) in the polytropic disk 
than in locally isothermal disks. S';>o have small but non- 
zero amplitudes near co-rotation, and their amplitude in 
the wave- like regions are at most ~ 20% of |S'o(ro)|. 

In the wave- like region, l^ool can be comparable or 
larger than l^o]. We found the solution in the wave re- 
gions more strongly affected by boundary conditions than 
in locally isothermal disks. 

We remark that for m — 5, 6, 5*0 no longer has 
largest disturbance amplitude around rp, because ra- 
dial confinement around co-rotation requires low m 
(jLin fc Papaloizoul l2011a[ ) , unless the vortensity mini- 
mum is very deep. At sufficiently large m (which de- 
pends on parameters), the modes are dominated by the 
wave-like region (much like the m = 5 mode in locally 
isothermal disks, see Fig. [3]). Boundary conditions are 
likely to play a role here, but they are not the vortex- 
forming RWI modes of interest. 
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Fig. 11. — Vertical velocity field for the m = 3 mode in the 
n = 1.5 fiducial polytropic disk. The slice is taken in the {r, z) 
plane at the azimuth = </<o. Arrows arc the perturbed velocity 
field projected onto this plane. 

6.2. Vertical structure 

We now examine the m = 3 mode in more detail. The 
flow in the (r, 0) plane is similar to the locally isother- 
mal disk. However, consistent with the previous section, 
vertical motion was found to be more prominent in the 
polytropic disk. 

As before we focus on the region r e [0.9, 1.1]. Fig. [TT] 
shows upwards vertical motion at the vortex core and is 
largest near z = H . The flow for z/ H < 0.5 and/or away 
from ro is essentially horizontal. The converging fiow 
pattern in Fig. [TT] is consistent with (ro,(?!)o) being an 
over-density. At the vortex core, upwards motion makes 
sense since the midplane is refiecting. It also implies an 
increase in disk thickness at (rg, 0o)- 

The background polytropic disk becomes thicker at rg 
(i.e. H varies on a local scale). Fluid moving into the 
vortex core finds itself in a region of larger vertical extent. 
Upwards motion enhances the disk thickness, consistent 
with enhanced pressure and with the RWI vortices being 
over-pressure regions. 

In the locally isothermal disk, it is difficult to di- 
rectly associate vertical motion with enhanced pressure 
as above, since the scale-height is prescribed to vary on 
a global scale and it remains unperturbed. Hence, verti- 
cal motion at (rg, 0g) was not seen in locally isothermal 
disks. 

We have also examined the vertical fiow in the poly- 
tropic disk for other m (< 4), but found similar fiow 
structure. This is unlike the locally isothermal disk which 
can display a range of vertical flow pattern depending on 
m. (Fig. [4]— [5]). This hints that there is a physical rea- 
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Fig. 12. — Vertical velocity field for the m = 3 mode in the 
n = 1.5 fiducial polytropic disk. The slice is taken in the {(f), z) 
plane at radius r = ro- Arrows are the perturbed velocity field 
projected onto this plane. 
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Fig. 14. — Same as Fig. 1131 but as a function of bump amplitude 
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Fig. 13. — Effect of h on three-dimensionality of the co-rotation 
flow (solid). Also shown are the normalized amplitude of 52 in 
this region (dashed) and the growth rates (dotted). The increase 
in growth rates with h is exp ected because the RWI is driven by 
pressure forces IILi et al.ll200CII ). 

son why polytropic disks tend to have positive vertical 
velocity at tq. We return to this point later. 

Lastly, Fig. [T^ shows the vertical flow in the (0, z) 
plane at r = tq. The flow is similar to that in the locally 
isothermal disk (Fig. ^ except that the region (f> ^ (f>Q is 
not in vertical hydrostatic equilibrium. 

6.3. Effect of h and A 

We measure the three-dimensionality of the flow in the 
co-rotation region in the same way as in H5.41 but here 
the averages are taken over the finite vertical extent of 
the disk. 

Fig. [TBI— fin show results from calculations with vari- 
able h (at fixed A — 1.4) and variable A (at fixed 
h = 0.14), respectively. The range of growth rates are 
similar to the cases examined for the locally isothermal 
disk (see Fig. [7]). I7I and (|-'^2|) also behave similarly. 

As in locally isothermal disks, Fig. [T31— [l4l shows that 
the three-dimensionality of the flow decreases with in- 
stability strength, but less rapidly in polytropic disks. 
Overall, (|i5v2|)/(|(5i'r|) does not vary much, consistent 
with our findings that the vertical flow structure, such 
as Fig. [n] — [T2I to be generic. Such plots are qualita- 
tively similar across the range of h and A considered. 
The vertical flow at the vortex core is always upwards. 

When the spatial average is taken over r £ [0.98, 1.02], 
we found (Sv^) /{\5vr\) maximizes at ft, = 0.16 for fixed 
A and a.t A — 1.6 for fixed h. However, its val- 
ues are of similar size: {6vz) / {\Svr\) — 0.44 — 0.65 and 
{Svz) / {\Svr\) — 0.53 — 0.65 for variable h and A, respec- 



tively. A reason for such insensitivity is that the above 
calculations have fixed polytropic index n, thereby fix- 
ing the fiuid properties. Below, we show that varying n 
affects the vertical flow. 

6.4. Other polytropic indices 

The polytropic index n not only affects the magnitude 
of the bump in the background disk thickness but also the 
compressibility of the fiuid. An isothermal fiuid can be 
considered a polytrope with n ^- 00 and is highly com- 
pressible, while n — corresponds to an incompressible 
fluid. Thus increasing n also increases compressibility. 

For polytropic disks we identified vertical flow at the 
vortex core. Here, we focus on this region and take radial 
averages over r e [0.98, 1.02]. Fig. [TSlshow calculations 
for n G [1,2.4]. As n decreases, instability strength in- 
creases and the vertical flow at ro noticeably increases, 
so the motion becomes more three-dimensional. This is 
qualitatively different from varying h or A, where the 
vertical fiow at the vortex core remain of similar size. 

At the CO- rotation radius, which is close to ro, the ver- 
tical velocity is 



\Svz 






{W 



1 



ZS2 



iHi 



(46) 



Hq is constant for fixed h. The factor |(4ri^ — 1)/7| de- 
creases with decreasing n, which by itself would reduce 
the vertical velocity. Fig. [15] shows this is not the case. 
The increase in 15*21 with decreasing n overcomes this 
effect. 

In Fig. [injwe compare the flow in the (r, z) plane be- 
tween n = 1 and n = 2. As previously remarked, the 
flows share the same qualitative feature: converging to- 
wards ro with upwards motion near ro. However, for 
smaller n (stronger instability), upwards motion is con- 
centrated at ro whereas for larger n (weaker instability) 
there is also upwards motion away from the vortex core. 
The latter was also seen for locally isothermal disks, con- 
sistent with larger n being more compressible. 

A larger vertical velocity at ro with decreasing n is 
consistent with variable compressibility. First note that 
IS'21 ^ l^oj in r e [0.9,1.1] so the perturbed enthalpy, 
radial and azimuthal velocities are all dominated by 
5*0, which gives converging flow towards the vortex core 
where there is enhanced pressure or density. We may 
then ask what vertical motion at ro is compatible with 
this 2D perturbed flow, as implied by 5*0? 
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Fig. 15. — Dependence of the vertical flow at the vortex core on 
the polytropic index n (soHd). The bump amphtude is fixed to 
A. = 1.4 and h = 0.14. Also shown are the growth rates (dotted) 
and amplitude of 52 (dashed). The mode is m, = 3. 
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Fig. 16. — Comparison between vertical velocity (contour) in a 
disk with polytropic index n = 1 (top) and n = 2 (bottom). Arrows 
indicate the velocity field projected onto this plane. 

At the vortex core (ro,(?!)o), the Hnearized continuity 
equation is approximately 

dt{Sp/p) -- -V • Sv - Svzdzlnp, 

where the S quantities are regarded as real. If the fluid 
is highly compressible (large n), then the density at the 
vortex core may increase with vertical motion playing no 
role. That is, the divergence term on the RHS dominates 
over the second (V • 5v itself dominated by horizontal 
velocities). 

However, if the fluid is made less compressible (de- 
creasing n), so that V • (5?; is reduced in magnitude, 
then the fluid at (rgjijio) should move upwards so that 
—Svzdz Inp > contributes to increasing the density. For 
n ^ 1, the fluid becomes incompressible so that V • 5v 
is negligible. Then the density can only increase by the 
fluid moving upwards, increasing the disk thickness and 
accommodating more material. 

It is important to note that in the above argument, 
we deduced vertical motion by imposing the 2D solu- 



tion in the three-dimensional disk. Effectively, we re- 
garded S'o is a source for 5*2, and that 5*2 has no back- 
reaction on S'o. This interpretation may not work for 
general disturbances, however. Here it is justified by 
the fact that |iS'2| ^ l^ol from the numerical calcula- 
tions. Calculations where the disk is truncated by set- 
ting ri — 0.7, To — 1.3, thereby excluding the wave-like 
regions in Si, show similar upwards motion. This indi- 
cates that S'o induces S2 locally. 

7. DISKS WITH k2 < 

IMeheut et al.l (|201 (y) performed the first nonlinear hy- 
drodynamic simulations that showed evidence for the 
RWI in a 3D polytropic disk. Their fiducial calculation 
showed the development of a m = 1 anti-cyclonic vortex 
which survived many orbits. 

Indeed, the consideration of polytropic disks in this 
paper was originally inspired by these simulations, but it 
turns out that the disk model employed bv lMeheut et ahl 
(|2010f ) has a region where k ^ < 0. Motivated b y this 
feature, in this section we use IMeheut et al.l ()201Clf )'s disk 
model to explore the 3D RWI when Kg < 0. We find 
that the RWI can be quite different to those described 
previously (where k^ > everywhere). 

It is straight forwar d to adapt our setups to models 
used bv lMeheut et al.l . They considered a n = 1.5 poly- 
tropic disk, occupying r G [^,^0] = [1,6], and specified 
the midplane density to be a power law {po oc r~^''^) with 
a Gaussian bump. Their bump in midplane density has 
the same functional form as that used for surface density 
in our models (Eq. [5]), so ^ is now the bump amplitude 
in midplane density. The bump is located at ro = 3 with 
width A?" = O.lri. The calculations presented below em- 
ployed Nr = 768 grid points, on account of the larger 
disk compared to previous models. 

We will consider the m = 1 mode below. Calcula- 
tions were done for m < 6, which gave similar growth 
rates when Kq < 0, but provided A is chosen to ensure 
Kq > 0, then higher m modes become dominant (e.g., 
with A = 1.15, TO = 5 had the highest growth rate). The 
latter is qualitat ively consistent with very recent numer- 
ical simulations (jMeheut et al.ll2012al . see also Appendix 

mi). 

When Kq > 0, we find similar flow structure to that 
described previously. Having applied the linear calcu- 
lations to a different disk model and recovering similar 
results gives us confidence in the robustness of the RWI 
to develop 3D. 

7.1. TO = 1 modes 



In their fiducial setup, IMeheut et al.l adopted a bump 
amplitude of ^ = 1.4. This results in k^ = at 
r ~ 0.99ro, l.Olro and Ko < (which is also reflected 
in their Fig. 9). The disk is therefore unstable to local 
axisymmetric perturbations (jChandrasekhadflQGlh . 

Interestingly, for A = lA we found a, m = I mode with 
large growth rate, I7I = 0.36iloi almost twice the largest 
growth rates found previously. Below, we examine this 
solution along with a case with A = 1.3, which has k^ > 
everywhere and growth rate I7I = 0.0517(j3. 

Despite A being similar, the m = 1 growth rate for 
A = 1.3 is much smaller than that for A = 1.4. For A — 

^ This is comparable to the nonlinear simulation. 
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Fig. 17. — Radial functions Si for the polytropic disk model 
with a bump in midplane density of amplitude A = lA (top, 
IMeheut et aLF s fiducial setup) and with a bump amplitude of 
A = 1.3 (bottom). For A = 1.4, ftg < 0. 

1.3 we did not find other m = I modes with growth rates 
similar to the m = 1 mode in A — 1.4. Furthermore, for 



A= lA the quantity D 
ro- 



'? — a^ almost vanishes near 



^\n{\D\IVL\) =4x 10-^ 

which occurs at r = 1.002ro. For A = 1.3, the value 
above is 0.14. 

Fig. [T71 compares the 5; functions for the cases above. 
While the double-peak in Sq for ^ = 1.3 w a s also found 
in previous sections and also by iLi et al.l (|200Cl[ ). it is 
absent in A ~ 1.4. The dominant 3D mode is ^2, but 
it is significantly larger in ^ = 1.4 than in ^ = 1.3. 
This indicates the vertical flow will also be qualitatively 
different. 

Fig. [T8l shows the flow pattern at = 0o for A ^ 1.4. 
This result is very different from that for A = 1.3, which 
share the same features as our previous setups with k^ > 
(e.g. Fig. [TT|) . Note that while the Si behave smoothly 
across rg (Fig. 1171) . numerical evaluation of Svr involves 
a division by D, which is very small near vq for A — 1.4. 
Thus, horizontal velocities may be subject to numerical 
artifacts at vq. Despite this, the direction of radial flow, 
being inwards for r < tq and outwards fo r r > ro with 
a shar p transition at rg, was also found in IMeheut et al.l 
(|2OT0l their Fig. 11). 

Neither A produced vertical flow consistent with that 
in IMeheut et al.l ()2010D where strong downwards flow at 
tq were identified with rolls excited on either side. By 
contrast the linear solutions have upwards motion and 
there is no vortical motion in the (r, z) plane. 

Despite using the same disk models, several factors 
may have contributed to th e discrepancy between the 
linear calculation above and IMeheut et al.f s simulation. 
These include the treatment of the vertical domain, non- 




FlG. 18. — Verti cal velocity field for a m = 1 mode in 
IMeheut et al.l l|201CI 1's disk model with A = 1.4, which results in 
in K^ < at ro and D — > there as well. The arrows are the 
perturbed velocity field projected onto this plane. This result is 
qualitatively different to modes with k'^ > (see figures in Q. 

linearities (H. Meheut, private communication) and in- 
teraction with other m modes in the simulation which 
cannot be treated in linear theory. 

There may also be numerical issues in our linear cal- 
culation because of k^ < 0. The RWI is associated with 
the term oc {l/a)dr{pofl/D) and its disturbance is local- 
ized about rg. This term is ex 1/k^, which diverges when 
K^ — >• near rg. Also because Kg < 0, it allows D ^ 
at co-rotation as well. We have performed calculations 
with lower spatial resolution, so that numerically k^ and 
D have larger deviations from zero, but we found simi- 
lar eigenfrequcncies and flow patterns to the case shown 
above. Wc will further comment on RWI modes with 



8. SUMMARY AND DISCUSSION 

In this paper, we have examined the linear stability of 
three-dimensional, vertically stratified and radially struc- 
tured disks. Our cal c ulatio ns are 3D analogs of those 
presented by ILi et "aTl (|200Clf ). in which the Rossby wave 
instability was studied in razor-thin disks. In order to 
simplify the problem, we assumed the perturbed hydro- 
dynamic quantities have vertical dependence that can be 
decomposed into Hermite or Gegenbauer polynomials. 
Our conclusions therefore apply to such perturbations 
only 

Our numerical calculations confirm the RWI persists 
in 3D. For ease of discussion below, we denote the full 
linear solution schematically as 

X = Y{r) + AY{r,z), 

where Y is the z-independent part of the solution and 
Ay is the part that also depends on z. 

8.1. Validity of 2D 

We find the RWI growth rate I7I can be accurately de- 
termined from the 2D problem alone. In other words, in- 
stability is associated with Y{r). In the region of interest 
— the vortensity minimum — where vortex-formation is 
expected, we find |Ay| <^ \Y\ so that enthalpy, radial 
velocity and azimuthal velocity perturbations have es- 
sentially no z-dependence. 

In fact, weak z-dependence is expected from ear- 
lier studies of accretion tori. For slender tori, 
iPapaloizou fc Pringld ()1985D demonstrated the existence 
of low m unstable modes with weak z-dependence. 
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iGoldreich et al.l ()1986l ) also justified the use of height- 
averaged equations for calculating modes in narrow tori, 
for which vertical hydrostatic equilibrium was assumed. 
Although we considered radially extended disks, their re- 
sults should apply here because the low m RWI modes, 
relevant to vortex-formation, have largest disturbance as- 
sociated with a narrow r e gion about the density bump. 
More recently, lUmurhaiil ()201Cl[ ) reproduced the RWI in 
approximate three-dimensional disk models, in which 
horizontal velocities have no vertical dependence. Our 
numerical results are therefore supported by analytic 
studies above. 

The 2D solution, Y, imply anti-cyclonic motion as- 
sociated with over-densities, thus we expect the RWI 
will lead to columnar vortices in 3D. The survival of 
vortices in 3D is then an imp ortant issue because they 
may be subject to instabilities (jLesur fc Papaloizoul2009L 
l2010f l. On the other hand, if there is a continuous source 
of vortensity extremum, such as disk-planet interaction, 
then vortex-formation via the RWI could be maintained. 

8.2. Vertical motion 

Although I AY I is small in the co-rotation region, it 
is nevertheless non-zero. This implies vertical motion 
growing on dynamical time-scales, making the flow in 
the co-rotation region three-dimensional. We found the 
nature of the vertical flow is affected by the equation of 
state. 

In polytropic disks the vortex core (ro,0o) always in- 
volve upwards motion. For fixed polytropic index n, 
there is limited variation in the magnitude of vertical 
flow with respect to instability strength. However, if the 
fluid is made less compressible by lowering n, then verti- 
cal motion at the vortex core increases. 

This result motivates us to interpret vertical motion 
arou nd co-rotation as a pe rturbation to the 2D solu- 
tion (jGoldreich et al.l I1986D . Recall that Y is the so- 
lution to the vertically integrated system. It signifies 
non-axisymmetric enhancements in surface density at the 
bump radius. This characteristic feature is unchanged by 
the addition of AY to the 2D solution. We then ask how 
should the disk respond in the vertical direction. 

The polytropic disk thickness is directly related to 
the surface density (Eq. [12]). Enhancement of the sur- 
face density therefore imply enhancement in disk thick- 
ness, so fluid at (ro,(/)o) moves upwards. If we look in 
the ((/>, z) plane at tq, the disk thickness becomes non- 
axisymmetric. T his has already been observed in nonlin- 
ear simulations (jMeheut et al.l [2011 al ibi). In these newer 
simulations, the authors indeed find upwards motion in 
anti-cyclonic vortices. 

Note that the polytropic disk thickness becomes less 
sensitive to surface density as n is increased (Eq. [T^j). 
For n -> oo the disk thickness is independent of surface 
density and there is no need for fluid to move vertically 
in order to achieve a surface density increase. In this case 
there is no preference for vertical velocity of a particular 
sign at (ro,(/)o). Since the fluid behaves isothermally as 
n — ^ oo, the above is consistent with our observation 
that locally isothermal disks have little vertical motion 
right at the vortex core. In Appendix IB. 2! we consider 
a polytropic disk calculation with n = 8 to check for 
consistency. 



8.3. RWI with k2 < 

We briefly examined the linear 3D RWI in disk 
where k^ becomes negative at the density bump. This 
was inspired by the 3D RWI simulations presented in 
iMeheut et ail (|2010[ ). where the disk model had k^ < 0. 
In this setup we found a to = 1 linear mode with large 
growth rate and qualitatively different to modes in disks 
with K^ > everywhere. In neith er case d id we reproduce 
the vertical flow seen in Mehe ut et al.l ([2010), namely 
downwards flow at the vortex center. 

Most discussions of non-axisymmetric disk instabil- 
ities have assumed k^ > everywhere, including 
iLovelace et al.l (|1999 ^ 's original description of the RWI, 
so that Rayleigh's criterion for stability against local ax- 
isymmetric perturbations is satisfled. 

The RWI has been shown to exist for Kg < but its 
properties ap pear diffe r ent to those in disks with Kg > 0. 
For example, iLi et al.l ()2000l )'s linear calculations indi- 
cate a non-smooth change in growth rate as Kg becomes 
negative (the 'HG B' case in their Fig. 9). In nonlinear 
2D simulations bv .Li et al.l (|2001l ). the RWI also evolves 
differently depending on whether the growth rate is low 
(|7| - O.lf^o and k^ > 0) or high (I7I - 0.317o and 
Kq < 0). Note the latter case has I7I close to that found 
in our calculation. We therefore expect the RWI to dif- 
fer in 3D depending on sgn(KQ). This is apparent by 
comparing our results with K q > t o those with Kg < 0. 

Thus, while I Meheut et all (|2010f) is the flrst demon- 
stration of the 3D RWI, it should be kept in mind that the 
disk model has Kg < 0. An understanding of such modes 
in 3D is of theoretical interest, but it is unclear whether 
or not protoplanetary disks will develop sufficiently large 
press ure gradients to render k^ < ([Yang fc Menoul 
[2OT0h . 

8.4. Outstanding issues 

The main goal of our study is to demonstrate the linear 
RWI in 3D and to identify the nature of associated three- 
dimensional flow structure around co-rotation. However, 
our study is subject to several caveats which should be 
clarified in future work. 

8.4.1. Baroclinic effects 

One issue is that our locally isothermal basic states are 
not in true equilibrium, because we approximated the ro- 
tation profile to be z-indcpendcnt (Eq. [5|). Initializing a 
full hydrodynamic simulation this way might boost radial 
velocities because of the inexact radial momentum bal- 
ance. In order for the angular velocity to be strictly in- 
dependent of z, we must set H ex r^^^, whi ch is the glob- 
ally isothermal disk already considered bv M eheut et al.l 
fj2012b.) . We do not expect this to make a difference from 
our disks with H (x r, because the RWI is driven by local 
variations in disk structure and its disturbance is radially 
confined. We check this in Appendix IB. 21 

Another justification is that for a thin, smooth disk 
(^ = 1) with H = hr, the angular velocity is 



Jl(r, z) = Clk 



h^ 
2 



2i/2 ; 



(47) 



([Tanaka et al.l 120021 ) . The difference in angular speed 
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between the gas at the niidplane and gas at z is then 

Ar! = |f](r,z)-0(r, 0)1 =/.2('_f_') rife (48) 

(a radial density bump does not contribute to this dif- 
ference). Since the gas is contained within a few scale- 
heights, we have Ail/fife = 0{h^). Because /i <C 1, verti- 
cal shear should be unimportant if the dynamics of inter- 
est operate on much faster time-scales, as can be the case 
for the RWI with growth rates ~ hftk- That is, the vor- 
tical perturbation grows much faster than it is sheared 
apart by Afi. We have begun preliminary nonlinear sim- 
ulations which confirms vortex formation via the RWI in 
a locally isothermal 3D disk with constant aspect-ratio 
( Lin 2012, in prepara t ion). 

iKnobloch fc Spruitl ([1986^ have pointed out the possi- 
bility of baroclinic instability in the case of fi = n,{r, z), 
when there are radial variations in temperature on the 
scale of local scale-heights. This condition is not met 
in our locally isothermal disk models because the sound- 
speed varies on a global scale. In more realistic disk mod- 
els, one might expect that a density bump also involves 
local temperature variations. Baroclinic effects may then 
become important. On the other hand, the RWI may 
also be enhanc ed because of local temperature gradients 
(jLi et al.ll200C0 . Having fl = fl{r,z) means solving the 
linearized equations as a PDE eigenvalue problem, which 
is not simple. 



ally takes place. Distant radial boundaries do not af- 
fect the dynamics in this region significantly (as checked 
numerically). However, it is clear that far away from 
co-rotation, three-dimensional effects become increas- 
ingly important. This is seen in the polytropic disk as 
|Ar| -- |r| towards the disk boundaries (Fig. ^. Dis- 
turbances associated with the RWI are therefore three- 
dimensional beyond the Lindblad resonances. In order 
to study these regions, more physically realistic radial 
boundary conditions are needed. 

Around co-rotation the RWI is a global disturbance in 
z, so the upper disk boundary conditions could be im- 
portant. The use of orthogonal polynomials means we 
simply impose a regularity condition at the upper disk 
boundary (Q. This method of solution does not allow 
us to explore the effect of other vertical boundary con- 
ditions. Again, such a study involves a PDE eigenvalue 
problem, but can reveal to what extent the dominance of 
the 2D solutions found here are influenced by the specific 
decompositions employed. This will be the subject of a 
follow up paper. 

Nevertheless, we can make some speculations based on 
results here. The vanishing density at the polytropic disk 
surface is likely to provide a reflective upper boundary. 
This effect may be important. It might reduce the growth 
of the RWI if it remains predominantly a 2D disturbance, 
because the 2D solution alters the surface density, which 
is directly related to the disk thickness for a polytrope, 
but the disk thickness cannot change. 



8.4.2. Boundary effects 

We have restricted our attention to the co-rotation 
region because this is where vortex-formation eventu- 



I thank H. Meheut for useful discussions and clarifica- 
tion of their simulation results. I also thank O. Umurhan 
for comments on the first version of this paper. 
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APPENDIX 

EXPLICIT EXPRESSIONS FOR THE LINEAR OPERATORS 

Locally isothermal disks with constant aspect-ratio 

For locally isothermal disks with H ^ hr with h being a constant and fl taken to be a function of radius only, the 
operators governing the linear problem are given by 



Bi = 
Ci = 



2mrfl d f cz 
a dr 

2 d \ 
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'^d-r- 



Amfll 



(Al) 
(A2) 
(A3) 



We have ex pressed the operators in terms of surface density S so the above may be seen to be equivalent to Eq. 21 in 
iTanaka eF al. (2002) when their parameter fi = dlnH/d\nr is set to unity. 

These equations are approximate because we ignored terms proportional to dzCl in the governing PDE from which 
Eq. lAll — [X3I are derived. These terms are non- vanishing for exact equilibrium if the sound-speed varies with radius, 
but for a thin disk dz^ oc h^ ^ 1 so we expect the m to be small. It is worth neglecting them in favor of the one- 
dimensional operators above, which are much simpler. ITanaka et al.l (|2002f ) gives a more general equation for the linear 
problem which includes dz^. Their Eq. 11 shows that Szil contributes to the coefficient of dzW as 



dW 
dz 



z m on 



dW 



m dz 



mh'^qflii 



(A4) 



where dzfl is evaluated using ITanaka et aD 's Eq. 4 and q = —d\ncs/d\nr — 0.5 for disks with constant aspect-ratio 
(equivalent to Eq. 05] in ^8.4.ip . Near co- rotation the magnitude of the second to first term is 



mh q 



rik 



mK^q 

\ll^k\ 



(A5) 



For the fiducial case va^^ m — i, h — 0.07 and |7/5^fc| — 0.057, this ratio is 0.13. We typically find |7/5^/£| = 
the second term is a factor mhq <C 1 smaller than the first for low m modes. Neglecting it (to arrive at Eq. 
is a self-consistent treatment. 
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Polytropic disks 
For polytropic disks, we find it most convenient to express the linear operators as 
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We have used the midplane density po, but it is straight forward to express the above in terms of S using the relation 
E = poH{r)In. The form of the operators above are appropriate for numerical computations in the range of polytropic 
indices considered in this paper (n > 1 or A > 0.5). Numerical issues may arise for smaller indices because of the 
(/ + A — 2)~^ factor in Bi. For example, if A = (n = 0.5) and / = 2 this factor diverges. However, for n — 0.5 it is 
more natural to use Chebyshev polynomials of the first kind (Tj) for expansion in z. We have performed calculations 
with n = 0.5 using Ti, and found similar results to those presented here. 

SUPPLEMENTARY CALCULATIONS 
Improved simulations 

During the finishing stages of this paper, IMeheut et al.l (|2012aD published new simulations of the 3D RWI with 
improved numerical resolution. This simulation developed a m = 5 mode with growth rate I7I = O.17f2oj with upwards 
motion at anti-cyclonic vortex centers and downwards motion at cyclonic vortex centers, which are consistent with our 
fiducial polytropic disks (Sj6]). 

We were able to find a rn = 5 linear mode provided the bump amplitude A in the midplane density was chosen to 
ensure k^ > 0. Using A = 1.7, we find I7I = O.ISJIq for m = 5. This mode is shown in Fig. [THl Note that Sq is still 
localized about rg, despite the higher m than those considered in our fiducial calculations (which gave more global 
disturbances). This is because here the vortensity minimum is deep, with inm{K'^/il.^) ~ 0.1, so even high m modes 
can b e localized . The v ertical flow at the vortex core is upwards, as found previously. 

Me heut et al.l ()2012a[ ) actually employed A — 2, giving Kq 
with similar growth rate as their simulation. As Kq is more negative in their new simulation than in IMeheut et al.l 
(|2010f) . one possibility is that an axisymmetric disturbance develops early on, rendering k^ > then the usual RWI 
follows. For A = 1.7 we find linear growth rates peak at m = 8 with I7I = 0.2irio! but this is only marginally larger 
than TO = 5. Differences in the linear and nonlinear calculations, such as the treatment of vertical boundaries, may 
then account for observation of m = 5 in the simulations. 
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Fig. 19. — A linear m = 5 mode found in lMeheut et al] l l2012al )'s n = 1.5 polytropic disk model, but with a smaller bump amplitude 
than their simulation. Left: radial eigenfunctions Si normalized by ISoC^o)!. Right: vertical flow structure. 




Fig. 20. — A m = 3 mode in the standard ra = 8 polytropic disk (growth rate I7I = 0.055f2o). Contours of the real vertical velocity 
perturbation are shown. Arrows are the perturbed velocity field projected onto this plane. This figure is similar to locally iso thermal disks 
in that there is very little vertical velocity near the vortex core, and is unlike polytropic disks with smaller n (e.g Fig. 1161 which shows 
significant upwards motion near r = ro). 




Fig. 21. — A linear mode in the locally isothermal disk with the same parameter values as the lMeheut et aLl l|2012b')'s globally isothermal 
disk (m = 4, h = O.l, other parameters are the same as our fiducial case in ijsjl. Contours of the real vertical velocity perturbation is shown. 
Arrows are the perturbed velocity field projected onto this plane. The growth rate I7I = 0.20r2o is similar to lMeheut et al.i . The vertical 
flow is also consistent with their Fig. 3d, namely the vertical velocity vanishes near r = tq. 

Consistency check 

We describe calculations to check the consistency between locally isothermal and polytropic disks and against the 
globally isothermal disk presented in Meheut et al. (2012b). 

Noting that an isothermal disk is a special case of a polytropic disk in the limit of large n, we performed a polytropic 
disk calculation with n = 8, yt = 2.0, /i = 0.2. Fig. EUlshow that in this case vertical motion is much smaller than the 
horizontal flow in the co-rotation region, compared to smaller values of n discussed in ^6.41 This is consistent with our 
typic al results for loc ally isothermal disks were the vertical velocity vanishes at the vortex core. 

Me heut et al.l (J201 21J') solved the linear problem for globally isothermal disks. Their basic state with fi = J7(r) 
satisfy exact radial momentum balance but adopting such a profile for locally isotherma l disks is only approximate. 
We have performed a locally isothermal calculation with the sa me parameters as Meheut et alJ (|2012b) . The result 
is shown in Fig. [5TJ It shares the same vertical flow implied by iMeheut et al.l ()2012bO 's Fig. 3d around a maximum 
in the (real) density perturbation: 5vz > near r = 1.1, 5vz < near r = 0.9 and Svz ~ at r = tq. This suggests 
that a node in the vertical velocity at the vortex core is a generic feature for linear RWI modes in isothermal disks. A 
global temperature profile does not affect the 3D RWI significantly. 



